//
// DistributionQaFigures.cxx
//
// Created by TB on 12/5/11.
// Copyright 2011 BNL. All rights reserved.
//

#include <TDirectory.h>
#include <TF1.h>
#include <TMath.h>

#include <boost/foreach.hpp>

#include "tppmc/DistributionQaFigures.h"
#include "tppmc/Distributions.h"
#include "tppmc/DssFragmentation.h"
#include "tppmc/Enumerations.h"
#include "tppmc/FragmentationFunction.h"
#include "tppmc/Particle.h"
#include "tppmc/PartonContent.h"
#include "tppmc/PartonDistribution.h"
#include "tppmc/RandomSpin.h"

namespace tppmc {
   
  DistributionQaFigures::DistributionQaFigures(TDirectory* dir)
    : fuquark("fuquark","u-quark fragmentation function (pion +)", 100, 0.05, 1.)
    , fubquark("fubquark","ubar-quark fragmentation function (pion +)", 100, 0.05, 1.) 
    , fdquark("fdquark","d-quark fragmentation function (pion +)", 100, 0.05, 1.) 
    , fdbquark("fdbquark","dbar-quark fragmentation function (pion +)", 100, 0.05, 1.) 
    , fsquark("fsquark","s-quark fragmentation function (pion +)", 100, 0.05, 1.) 
    , fsbquark("fsbquark","sbar-quark fragmentation function (pion +)", 100, 0.05, 1.) 
    , fgluon("fgluon","gluon fragmentation function (pion +)", 100, 0.05, 1.) 
   
    , fkuquark("fkuquark","u-quark fragmentation function (kaon +)", 100, 0.05, 1.) 
    , fkubquark("fkubquark","ubar-quark fragmentation function (kaon +)", 100, 0.05, 1.) 
    , fkdquark("fkdquark","d-quark fragmentation function (kaon +)", 100, 0.05, 1.) 
    , fkdbquark("fkdbquark","dbar-quark fragmentation function (kaon +)", 100, 0.05, 1.) 
    , fksquark("fksquark","s-quark fragmentation function (kaon +)", 100, 0.05, 1.) 
    , fksbquark("fksbquark","sbar-quark fragmentation function (kaon +)", 100, 0.05, 1.) 
    , fkgluon("fkgluon","gluon fragmentation function (kaon +)", 100, 0.05, 1.) 
   
    , fpuquark("fpuquark","u-quark fragmentation function (proton)", 100, 0.05, 1.) 
    , fpubquark("fpubquark","ubar-quark fragmentation function (proton)", 100, 0.05, 1.) 
    , fpdquark("fpdquark","d-quark fragmentation function (proton)", 100, 0.05, 1.) 
    , fpdbquark("fpdbquark","dbar-quark fragmentation function (proton)", 100, 0.05, 1.) 
    , fpsquark("fpsquark","s-quark fragmentation function (proton)", 100, 0.05, 1.) 
    , fpsbquark("fpsbquark","sbar-quark fragmentation function (proton)", 100, 0.05, 1.) 
    , fpgluon("fpgluon","gluon fragmentation function (proton)", 100, 0.05, 1.) 
   
    , fuquarkc("fuquarkc","u-quark fragmentation function (pion +, Collins)", 100, 0.05, 1.,100,0.0,TMath::TwoPi()) 
    , fubquarkc("fubquarkc","ubar-quark fragmentation function (pion +, Collins)", 100, 0.05, 1.,100,0.0,TMath::TwoPi()) 
    , fdquarkc("fdquarkc","d-quark fragmentation function (pion +, Collins)", 100, 0.05, 1.,100,0.0,TMath::TwoPi()) 
    , fdbquarkc("fdbquarkc","dbar-quark fragmentation function (pion +, Collins)", 100, 0.05, 1.,100,0.0,TMath::TwoPi()) 
   
    , favff("favff","u-quark fragmentation function (pion +)", 100, 0.05, 1.) 
    , unfavff("unfavff","d-quark fragmentation function (pion +)", 100, 0.05, 1.) 
   
    , ratiofav("ratiofav","#Delta^{N}D_{fav}(z)/2D_{fav}(z)",80,0.1,0.9) 
    , ratiounfav("ratiounfav","-#Delta^{N}D_{unf}(z)/2D_{unf}(z)",80,0.1,0.9)
   
    , ratiofav2("ratiofav2","#Delta^{N}D_{fav}(z)/2D_{fav}(z)",100,0.0,1.0) 
    , ratiounfav2("ratiounfav2","-#Delta^{N}D_{unf}(z)/2D_{unf}(z)",100,0.0,1.0) 
   
    , transdist_upq("transdist_upq","x#Delta_{T}u(x)",100,0.0,1.0)
    , transdist_dnq("transdist_dnq","x#Delta_{T}u(x)",100,0.0,1.0)
    , transdist_sb("transdist_sb","x#Delta_{T}u(x)",100,0.0,1.0)
   
    , transversity_upq("transversity_upq","x#Delta_{T}u(x)",100,0.0,1.0)
    , transversity_dnq("transversity_dnq","x#Delta_{T}d(x)",100,0.0,1.0)
    , transversity_upb("transversity_upb","x#Delta_{T}ubar(x)",100,0.0,1.0)
    , transversity_dnb("transversity_dnb","x#Delta_{T}dbar(x)",100,0.0,1.0)
   
    , unpol_pdf("unpol_pdf", "", 100, -3, 0.)
   
    , sivers_upq("sivers_upq","x#Delta^{N}f^{(1)}(x)",100,-3.0,0.0)
    , sivers_dnq("sivers_dnq","x#Delta^{N}f^{(1)}(x)",100,-3.0,0.0)
    , sivers_upb("sivers_upb","x#Delta^{N}f^{(1)}(x)",100,-3.0,0.0)
    , sivers_dnb("sivers_dnb","x#Delta^{N}f^{(1)}(x)",100,-3.0,0.0)
    , sivers_sq("sivers_sq","x#Delta^{N}f^{(1)}(x)",100,-3.0,0.0)
    , sivers_sb("sivers_sb","x#Delta^{N}f^{(1)}(x)",100,-3.0,0.0)
    , sivers_glu("sivers_glu","x#Delta^{N}f^{(1)}(x)",100,-3.0,0.0)
   
    , sivers_upq_k("sivers_upq_k","x#Delta^{N}f^{(1)}(x,k)",100,0.0,1.1)
    , sivers_dnq_k("sivers_dnq_k","x#Delta^{N}f^{(1)}(x,k)",100,0.0,1.1)
    , sivers_upb_k("sivers_upb_k","x#Delta^{N}f^{(1)}(x,k)",100,0.0,1.1)
    , sivers_dnb_k("sivers_dnb_k","x#Delta^{N}f^{(1)}(x,k)",100,0.0,1.1)
    , sivers_sq_k("sivers_sq_k","x#Delta^{N}f^{(1)}(x,k)",100,0.0,1.1)
    , sivers_sb_k("sivers_sb_k","x#Delta^{N}f^{(1)}(x,k)",100,0.0,1.1)
    , sivers_glu_k("sivers_glu_k","x#Delta^{N}f^{(1)}(x,k)",100,0.0,1.1)
   
    , fffav2d("fffav2d","",100,0.,1.,100,-TMath::Pi(),TMath::Pi())
  {   
    Create(dir);
  }
   
  DistributionQaFigures::~DistributionQaFigures() {
  }
   
  void DistributionQaFigures::Create(TDirectory* dir) {
    mFigures.clear();
      
    TDirectory* initialDir = gDirectory;
    bool addDir = TH1::AddDirectoryStatus();
      
    TH1::AddDirectory(true);
      
    // Use the current directory if one was not specified.
    if(not dir and gDirectory)
      dir = gDirectory;
    if(dir) dir->cd();
      
    BOOST_FOREACH(TH1* h, mFigures) {
      h->SetDirectory(dir);
    } // BOOST_FOREACH
    //      gDirectory->ls();
      
    mFigures.push_back(&fuquark);
    mFigures.push_back(&fubquark);
    mFigures.push_back(&fdquark);
    mFigures.push_back(&fdbquark);
    mFigures.push_back(&fsquark);
    mFigures.push_back(&fsbquark);
    mFigures.push_back(&fgluon);
      
    mFigures.push_back(&fkuquark);
    mFigures.push_back(&fkubquark);
    mFigures.push_back(&fkdquark);
    mFigures.push_back(&fkdbquark);
    mFigures.push_back(&fksquark);
    mFigures.push_back(&fksbquark);
    mFigures.push_back(&fkgluon);
      
    mFigures.push_back(&fpuquark);
    mFigures.push_back(&fpubquark);
    mFigures.push_back(&fpdquark);
    mFigures.push_back(&fpdbquark);
    mFigures.push_back(&fpsquark);
    mFigures.push_back(&fpsbquark);
    mFigures.push_back(&fpgluon);
      
    mFigures.push_back(&fuquarkc);
    mFigures.push_back(&fubquarkc);
    mFigures.push_back(&fdquarkc);
    mFigures.push_back(&fdbquarkc);
      
    mFigures.push_back(&favff);
    mFigures.push_back(&unfavff);
      
    mFigures.push_back(&ratiofav);
    mFigures.push_back(&ratiounfav);
      
    mFigures.push_back(&ratiofav2);
    mFigures.push_back(&ratiounfav2);
      
    mFigures.push_back(&transdist_upq);
    mFigures.push_back(&transdist_dnq);
    mFigures.push_back(&transdist_sb);
      
    mFigures.push_back(&transversity_upq);
    mFigures.push_back(&transversity_dnq);
    mFigures.push_back(&transversity_upb);
    mFigures.push_back(&transversity_dnb);
      
    mFigures.push_back(&unpol_pdf);
      
    mFigures.push_back(&sivers_upq);
    mFigures.push_back(&sivers_dnq);
    mFigures.push_back(&sivers_upb);
    mFigures.push_back(&sivers_dnb);
    mFigures.push_back(&sivers_sq);
    mFigures.push_back(&sivers_sb);
    mFigures.push_back(&sivers_glu);
      
    mFigures.push_back(&sivers_upq_k);
    mFigures.push_back(&sivers_dnq_k);
    mFigures.push_back(&sivers_upb_k);
    mFigures.push_back(&sivers_dnb_k);
    mFigures.push_back(&sivers_sq_k);
    mFigures.push_back(&sivers_sb_k);
    mFigures.push_back(&sivers_glu_k);
      
    mFigures.push_back(&fffav2d);
      
    TH1::AddDirectory(addDir);
    gDirectory = initialDir;
    if(gDirectory) gDirectory->cd();
  }
   
   
  void DistributionQaFigures::Fill(unsigned ntimes,
				   FragmentationFunction& ff,
				   TH1& histogram,
				   const Parton& parton) {
    //      std::cout << "Filling " << histogram.GetName() << std::endl;
    for(unsigned i(0); i < ntimes; ++i) {
      std::vector<Hadron> hadrons = ff.Fragment(parton,
						::pow(2.5, 2));
      BOOST_FOREACH(Hadron& hadron, hadrons) {
	histogram.Fill(hadron.Z());
      }
    } // if
#if 0
    ff.Function()->SetParameter(1, 2.5);
    ff.Function()->SetParameter(2, parton.Type());
    ff.Function()->SetParameter(3, 0);
    ff.Function()->SetParameter(4, 1);
    ff.Function()->SetParameter(5, 1);
    histogram.FillRandom(ff.Function()->GetName(), ntimes);
#endif
  }
   
   
  //   frag->SetParameter(0, DSS);
  //   frag->SetParameter(1, ::sqrt(2.5)); // Scale Q
  //   frag->SetParameter(2, UQUARK);
  //   frag->SetParameter(3, LO);
  //   frag->SetParameter(4, 1.); // Species
  //   frag->SetParameter(5, 1.); // Charge
  void DistributionQaFigures::Fill(Distributions& d) {
      
    const int nFills(10000);
      
    FragmentationFunction& ff = *(d.Fragmentation());
    const int initialHadron(ff.HadronType());
    //      std::cout << ff.Name() << " initial forced hadron = " << initialHadron << std::endl;
    //      if(ff.Supports(PIPLUS)) {
    ff.SetHadron(PIPLUS);
    //         std::cout << "Plot for hadron " << ff.Hadron() << std::endl;
    Fill(nFills, ff, fgluon, Parton(0., 0., 1., 21));
    Fill(nFills, ff, fuquark, Parton(0., 0., 1., 2));
    Fill(nFills, ff, fubquark, Parton(0., 0., 1., -2));
    Fill(nFills, ff, fdquark, Parton(0., 0., 1., 1));
    Fill(nFills, ff, fdbquark, Parton(0., 0., 1., -1));
    Fill(nFills, ff, fsquark, Parton(0., 0., 1., 3));
    Fill(nFills, ff, fsbquark, Parton(0., 0., 1., -3));
    //      std::cout << "Plot for hadron " << ff.Hadron() << std::endl;
    // Favoured: up quark --> pi+
    Fill(nFills, ff, favff, Parton(0., 0., 0., 2));
    //      std::cout << "Plot for hadron " << ff.Hadron() << std::endl;
    // Unfavoured: down quark --> pi+
    Fill(nFills, ff, unfavff, Parton(0., 0., 0., 1));
    //      }
    //      BOOST_FOREACH(TH1* h, mFigures) {
    //         h->Scale(ff.Function()->Integral(0.05, 1.)/(nFills*((0.05-1.)/100.0)));
    //      }
      
    //      if(ff.Supports(KPLUS)) {
    ff.SetHadron(KPLUS);
    //         std::cout << "Plot for hadron " << ff.Hadron() << std::endl;
    Fill(nFills, ff, fkgluon, Parton(0., 0., 0., 21));
    Fill(nFills, ff, fkuquark, Parton(0., 0., 0., 2));
    Fill(nFills, ff, fkubquark, Parton(0., 0., 0., -2));
    Fill(nFills, ff, fkdquark, Parton(0., 0., 0., 1));
    Fill(nFills, ff, fkdbquark, Parton(0., 0., 0., -1));
    Fill(nFills, ff, fksquark, Parton(0., 0., 0., 3));
    Fill(nFills, ff, fksbquark, Parton(0., 0., 0., -3));
    //      }
    //      if(ff.Supports(PROTON)) {
    ff.SetHadron(PROTON);
    //         std::cout << "Plot for hadron " << ff.Hadron() << std::endl;
    Fill(nFills, ff, fpgluon, Parton(0., 0., 0., 21));
    Fill(nFills, ff, fpuquark, Parton(0., 0., 0., 2));
    Fill(nFills, ff, fpubquark, Parton(0., 0., 0., -2));
    Fill(nFills, ff, fpdquark, Parton(0., 0., 0., 1));
    Fill(nFills, ff, fpdbquark, Parton(0., 0., 0., -1));
    Fill(nFills, ff, fpsquark, Parton(0., 0., 0., 3));
    Fill(nFills, ff, fpsbquark, Parton(0., 0., 0., -3));
    //      }
    // plot the "favored" and "unfavored" pion fragmentation functions
      
    // Fill the transversity and Sivers distributions
      
    DssFragmentation dss;
      
    Parton upq(1., 1., 0.,  2);
    upq.cos_spin_depol = 1.;
      
    Hadron piplus;
    piplus.type = 211;
    piplus.px = 3.04863256751186806e-01;
    piplus.py = 0.;
    piplus.pz = 5.0;
      
    Hadron piminus(piplus);
    piminus.type = -211;
      
      
    for(int i(1); i <= ratiofav.GetNbinsX(); ++i) {
      piplus.z = ratiofav.GetBinCenter(i);
      double withCollins    = ff.Evaluate(upq, piplus, 2.4);
      double withoutCollins = dss.Evaluate(upq, piplus, 2.4);
      ratiofav.SetBinContent(i, fabs(withCollins - withoutCollins) / withoutCollins);
      //         std::cout << "with = " << withCollins << " without = " << withoutCollins << std::endl;
      //         std::cout << "z = " << piplus.z << " ratio = " << abs(withCollins - withoutCollins) / withoutCollins << std::endl;
      piminus.z = piplus.z;
      withCollins    = ff.Evaluate(upq, piminus, 2.4);
      withoutCollins = dss.Evaluate(upq, piminus, 2.4);
      ratiounfav.SetBinContent(i, fabs(withCollins - withoutCollins) / withoutCollins);
      //         std::cout << "with = " << withCollins << " without = " << withoutCollins << std::endl;
      //         std::cout << "z = " << piminus.z << " ratio = " << abs(withCollins - withoutCollins) / withoutCollins << std::endl;
    } // for
      
      
    PartonContent partons;
    RandomSpin spin(d.Unpolarised().get(), d.Transversity().get());
    for(int ti=0; ti<100; ti++){ 
      Parton p;
      p.x = (ti + 0.5) * 0.01;
      double q_squared = 2.4;
      d.Transversity()->Evaluate(p, q_squared, partons);
         
#if 0
      transversity_upq.Fill(p.X(), p.X() * spin.Generate(p, q_squared));
      transversity_dnq.Fill(p.X(), p.X() * spin.Generate(p, q_squared));
      transversity_upb.Fill(p.X(), p.X() * spin.Generate(p, q_squared));
      transversity_dnb.Fill(p.X(), p.X() * spin.Generate(p, q_squared));
#endif
      transversity_upq.Fill(p.X(), p.X() * partons[LHAPDF::UP]);
      transversity_dnq.Fill(p.X(), p.X() * partons[LHAPDF::DOWN]);
      transversity_upb.Fill(p.X(), p.X() * partons[LHAPDF::UBAR]);
      transversity_dnb.Fill(p.X(), p.X() * partons[LHAPDF::DBAR]);
      //         PARTONCONTENT transv;  
      //         TransversityDistribution(transv,  Params, parton_x1, q_squared);
      //         transdist_upq.Fill(parton_x1, parton_x1 * transv.up);
      //         transdist_dnq.Fill(parton_x1, parton_x1 * transv.down);
      //         transdist_sb.Fill(parton_x1, parton_x1 * transv.anti_strange);
    }
      
    for(int i(1); i <= fffav2d.GetNbinsX(); ++i) {
      for(int j(1); j <= fffav2d.GetNbinsY(); ++j) {
	double z   = fffav2d.GetXaxis()->GetBinCenter(i);
	if(z<ff.ZMin()) continue;
	piplus.z = z;
	double phi = fffav2d.GetYaxis()->GetBinCenter(j);
	TVector3 v(3.04863256751186806e-01, 0., 1.);
	v.SetPhi(phi);
	piplus.px = v.Px();
	piplus.py = v.Py();
	piplus.pz = v.Pz();
	double val = ff.Evaluate(upq, piplus, 2.4*2.4);
	fffav2d.SetBinContent(i, j, val);
      } // for
    } // for
      
    for(int ti=0; ti<100; ti++){ 
         
      // plot is in log(x)
      const double q_squared = 2.4; 
      const double l_parton_x1 = -3.0 + (ti+0.5)*(0.03);
      const double parton_x1 = pow(10.0,l_parton_x1); 
         
      Parton p;
      p.x = parton_x1;
         
      double sw(0.f);
         
      PartonContent unpol;
      d.Unpolarised()->Evaluate(p, q_squared, unpol);
         
      //         unpolarised(unpol, parton_x1, q_squared);
         
      unpol_pdf.Fill(l_parton_x1, parton_x1 * unpol[LHAPDF::UP]);
         
      p.py = 0.;
         
      sw = 0;
      for(int tj=0; tj<1000; tj++){
	float kperp = (tj+0.5)*(3.0/1000.0);
	p.px = kperp;
	PartonContent sivers;
	d.Sivers()->Evaluate(p, q_squared, sivers);
	sw += (sivers[LHAPDF::UP]
	       *(0.003)*kperp*kperp*(1.0/(4.0*0.938))*TMath::TwoPi())
	  *(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25); 
      }
      sivers_upq.Fill(l_parton_x1, sw * parton_x1);
         
      sw = 0;
      for(int tj=0; tj<1000; tj++){
	float kperp = (tj+0.5)*(3.0/1000.0); 
	p.px = kperp;
	PartonContent sivers;
	d.Sivers()->Evaluate(p, q_squared, sivers);
	sw += (sivers[LHAPDF::DOWN]
	       *(0.003)*kperp*kperp*(1.0/(4.0*0.938))*TMath::TwoPi()) 
	  *(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25); 
      }
      sivers_dnq.Fill(l_parton_x1, sw * parton_x1);
         
      sw = 0;
      for(int tj=0; tj<1000; tj++){
	float kperp = (tj+0.5)*(3.0/1000.0); 
	p.px = kperp;
	PartonContent sivers;
	d.Sivers()->Evaluate(p, q_squared, sivers);
	sw += (sivers[LHAPDF::STRANGE]
	       *(0.003)*kperp*kperp*(1.0/(4.0*0.938))*TMath::TwoPi()) 
	  *(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25); 
      }
      sivers_sq.Fill(l_parton_x1, sw * parton_x1);
         
      sw = 0;
      for(int tj=0; tj<1000; tj++){
	float kperp = (tj+0.5)*(3.0/1000.0); 
	p.px = kperp;
	PartonContent sivers;
	d.Sivers()->Evaluate(p, q_squared, sivers);
	sw += (sivers[LHAPDF::UBAR]
	       *(0.003)*kperp*kperp*(1.0/(4.0*0.938))*TMath::TwoPi()) 
	  *(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25); 
      }
      sivers_upb.Fill(l_parton_x1, sw * parton_x1);
         
      sw = 0;
      for(int tj=0; tj<1000; tj++){
	float kperp = (tj+0.5)*(3.0/1000.0); 
	p.px = kperp;
	PartonContent sivers;
	d.Sivers()->Evaluate(p, q_squared, sivers);
	sw += (sivers[LHAPDF::DBAR]
	       *(0.003)*kperp*kperp*(1.0/(4.0*0.938))*TMath::TwoPi())
	  *(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25); 
      }
      sivers_dnb.Fill(l_parton_x1, sw * parton_x1);
         
      sw = 0;
      for(int tj=0; tj<100; tj++){
	float kperp = (tj+0.5)*(3.0/1000.0); 
	p.px = kperp;
	PartonContent sivers;
	d.Sivers()->Evaluate(p, q_squared, sivers);
	sw += (sivers[LHAPDF::SBAR]
	       *(0.003)*kperp*kperp*(1.0/(4.0*0.938))*TMath::TwoPi())
	  *(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25); 
      }
      sivers_sb.Fill(l_parton_x1, sw * parton_x1);
         
      sw = 0;
      for(int tj=0; tj<1000; tj++){
	float kperp = (tj+0.5)*(3.0/1000.0); 
	p.px = kperp;
	PartonContent sivers;
	d.Sivers()->Evaluate(p, q_squared, sivers);
	sw += (sivers[LHAPDF::GLUON]
	       *(0.003)*kperp*kperp*(1.0/(4.0*0.938))*TMath::TwoPi())
	  *(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25); 
      }
      sivers_glu.Fill(l_parton_x1, sw * parton_x1);
    }
      
      
      
      
    for(int ti=0; ti<100; ti++){ 
         
      // plot is in kperp
         
      double q_squared = 2.4; 
      double parton_x1 = 0.1; 
      float kperp = (ti+0.5)*(1.1/100.0); 
         
      Parton p(kperp, 0., 0.);
      p.x = parton_x1;
         
      double sw(0.);
         
      PartonContent unpol; 
      d.Unpolarised()->Evaluate(p, q_squared, unpol);
         
      PartonContent sivers;
      d.Sivers()->Evaluate(p, q_squared, sivers);
         
      sw = sivers[LHAPDF::UP];
      sivers_upq_k.Fill(kperp, sw *parton_x1*(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25));	
         
      sw = sivers[LHAPDF::DOWN];
      sivers_dnq_k.Fill(kperp, sw *parton_x1*(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25));	
         
      sw = sivers[LHAPDF::UBAR];
      sivers_upb_k.Fill(kperp, sw *parton_x1*(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25));	
         
      sw = sivers[LHAPDF::DBAR];
      sivers_dnb_k.Fill(kperp, sw *parton_x1*(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25));	
         
      sw = sivers[LHAPDF::STRANGE];
      sivers_sq_k.Fill(kperp, sw *parton_x1*(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25));	
         
      sw = sivers[LHAPDF::SBAR];
      sivers_sb_k.Fill(kperp, sw *parton_x1*(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25));	
         
      sw = sivers[LHAPDF::GLUON];
      sivers_glu_k.Fill(kperp, sw *parton_x1*(1.0/(TMath::Pi()*0.25))*exp(-(kperp*kperp)/0.25));	
    }
      
      
      
      
    // Revert to initial value
    ff.SetHadron(initialHadron);
    //      std::cout << ff.Name() << " forced hadron reverted to " << ff.Hadron() << std::endl;
  }
   
  void DistributionQaFigures::Write() {
    BOOST_FOREACH(TH1* o, mFigures) {
      o->Write();
    }
  }
   
} // namespace tppmc
